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ABSTRACT 

Ground and space-based sky surveys enable powerful cosmological probes based on measurements of galaxy 
properties and the distribution of galaxies in the Universe. These probes include weak lensing, baryon acous¬ 
tic oscillations, abundance of galaxy clusters, and redshift space distortions; they are essential to improving our 
knowledge of the nature of dark energy. On the theory and modeling front, large-scale simulations of cosmic 
structure formation play an important role in interpreting the observations and in the challenging task of extract¬ 
ing cosmological physics at the needed precision. These simulations must cover a parameter range beyond the 
standard six cosmological parameters and need to be run at high mass and force resolution. One key simulation- 
based task is the generation of accurate theoretical predictions for observables, via the method of emulation. Using 
a new sampling technique, we explore an 8-dimensional parameter space including massive neutrinos and a vari¬ 
able dark energy equation of state. We construct trial emulators using two surrogate models (the linear power 
spectrum and an approximate halo mass function). The new sampling method allows us to build precision emula¬ 
tors from just 26 cosmological models and to increase the emulator accuracy by adding new sets of simulations in 
a prescribed way. This allows emulator fidelity to be systematically improved as new observational data becomes 
available and higher accuracy is required. Finally, using one ACDM cosmology as an example, we study the de¬ 
mands imposed on a simulation campaign to achieve the required statistics and accuracy when building emulators 
for dark energy investigations. 

Subject headings: methods: statistical — cosmology: large-scale structure of the universe 


1. INTRODUCTION 

Cosmology has developed rapidly in the last few decades. 
The more qualitative original picture of the evolution of the 
Universe has evolved into a well-tested six-parameter “Stan¬ 
dard Model of Cosmology”. Remarkably, this simple model de¬ 
scribes all current observations at a level of accuracy of roughly 
a few percent, cementing a consistent picture of the Universe 
from a diverse set of cosmological probes. While this progress 
has truly revolutionized our knowledge of the Universe, a real 
understanding of its two main constituents, dark matter and 
dark energy, remains elusive. 

To make progress towards solving these mysteries, as well 
as shedding light on other questions such as the sum of neu¬ 
trino masses and the nature of primordial fluctuations, new 
ground- and space-based missions have been carried out (so- 
called Stage II experiments), are ongoing (Stage III), and are 
under construction or planned (Stage IV). These include the 


Baryon Oscillation Spectroscopic Survey (BOSS) (Schlegel, 
White, & Eisenstein| 2009 !, WiggleZ ( |Drinkwater et al.|2010| , 
the Dark Energy Survey (DE5 'll, the Large Synoptic Survey 


Telescope (LSST) (Abell et al. 


troscopic Instrument (DESP i, the Wide Field Infrared Survey 


2009} , the Dark Energy Spec- 


Telescope (WFIRST) (Spergel |2015 1 , and Euclid (Refregier et 


1 https://www.darkenergysurvey.org/ 


2 http://desi.lbl.gov/ 


al.|2010[ >. The major science targets of these missions specific 
to dark energy are measurements of the large scale structure 
in the Universe on different scales (along with supernova dis¬ 
tance measurements). Weak lensing, baryon acoustic oscilla¬ 
tions (BAO), abundance of galaxy clusters, and redshift space 
distortions are the most prominent probes aimed at investigat¬ 
ing the nature of dark energy. The upcoming surveys promise 
to obtain measurements of the main cosmology parameters at 
the 1% level accuracy if systematic errors can be sufficiently 
controlled. 

To reach the desired accuracies, the physics of matter clus¬ 
tering on significantly nonlinear scales must be accurately mod¬ 
eled. Currently, the only way to do this is via expensive numer¬ 
ical simulations. Depending on the scales of interest and the 
required accuracy, different physical effects must be accounted 
for. On large scales relevant to, e.g., BAO, gravity-only sim¬ 
ulations are likely sufficient for interpreting ongoing and up¬ 
coming measurements. For weak lensing, current surveys can 
be modeled with pure gravity simulations (see, e.g., the discus¬ 
sion in Heitmann et al. (20101) but the interpretation of future 


surveys will require a treatment of baryonic effects (e.g., gas 
dynamics, feedback) and associated subgrid modeling. 

Besides adding hydrodynamics effects and other modeling 
ingredients to the simulations, the cosmological parameter space 


has to be simultaneously widened. As argued in Heitmann 
et al. (2010|>, current observations do not have enough infor- 
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mation to constrain cosmological parameters beyond wCD.M 
scenarios - a mild extension of ACDM - at high accuracy. 
In the near future, however, it will be important to go be¬ 
yond wCDM by including, e.g., dynamical dark energy mod¬ 
els and at the same time accounting for “secondary” effects, 
currently too small to reliably extract from the measurements 
beyond setting upper or lower limits (such as a lower limit 
for m v , the neutrino mass sum). The work here is aimed at a 
research program that covers the following eight-dimensional 
space: pg = {oj m ,uJb,crs,h,n s ,Wo,w a ,u] v }. (Details are found in 
Section^) 


1.1. Simulation Campaigns 

If only a small number of simulations are required, the size of 
the parameter space is essentially irrelevant. This is, however, 
not the case in cosmology. Eventually one wishes to infer the 
values of the parameters from a set of observational measure¬ 
ments. Commonly used procedures for doing this all fall under 
the class of statistical inverse problems, where a large number 
of evaluations of the forward model (i.e., results from many 
simulations, at different specified parameter values) are neces¬ 
sary in order to use Markov chain Monte Carlo (MCMC) or 
other related methods. This may involve hundreds of thousands 
or even millions of forward model evaluations. (A similar issue 
is encountered in covariance estimation, where a large number 
of simulations are also needed.) 

At first sight this problem may appear to be computationally 
overwhelming but it is in fact possible to cover large parameter 
spaces in an efficient manner with a drastically reduced num¬ 
ber of simulations, provided certain smoothness conditions are 
met. The last requirement is largely satisfied in cosmological 
inverse problems, and this has allowed us to develop the “Cos¬ 
mic Calibration Framework” ([Heitmann et al.||2006| Habib et 


aL|2007[ Heitmann et al.|2009[ Lawrence et ai.||2010 Higdon 


et al.pOlO Heitmann et al.||2014 1 . The major components of 

the framework are: 


Sophisticated sampling methods that provide optimal 
coverage of the parameter space. In previous work we 
used a symmetric Latin hypercube design, as explained 
in Heitmann et al. ( 2009} (for a more technical presen¬ 
tation, see Santner et al.||2003) . The new strategy pre¬ 
sented here allows us to begin with a small set of cos¬ 
mological models from which to construct a first set of 
emulators. These simulations are then systematically 
augmented to improve emulation accuracy. This defines 
a notion of strong convergence (Cf. Section [3} as op¬ 
posed to the weak convergence previously demonstrated 
in |Habib et all ( |2007} . 


• High-dimensional interpolater based on Gaussian pro¬ 
cess modeling to build a prediction scheme (the emula¬ 
tor) for a diverse set of cosmological statistics of inter¬ 
est, e.g., the fluctuation power spectrum, the halo mass 
function, lensing statistics, etc. 


• Emulator-based sensitivity analysis, characterizing the 
dependence of the cosmological statistics on different 
cosmological and modeling parameters. 


• Final calibration step, which combines observational 
data and theoretical predictions to extract cosmological 
information from the data via MCMC methods. 


We will focus on the first three steps and construct precision 
emulators for the power spectrum and the mass function; we 
also discuss extensions to other already existing emulators, e.g., 
for the concentration-mass relation (|Kwan et al.|2013a| and the 


galaxy power spectrum and correlation function (Kwan et al. 
1 2017 5) below. 

In this paper, we narrow our focus to work that directly im¬ 
pacts gravity-only simulations. The purpose here is to demon¬ 
strate the success of the basic methodology. Over time, we 
expect to include direct modeling of baryonic effects as well 
as indirect approaches based on post-processing gravity-only 
simulations. As has been mentioned previously, and discussed 
in several papers (White 2004| Zhan & Knox|2004 |Jing et ah 


2006] Rudd et al. 2008; van Daalen et al.||2011} , deriving re 

liable predictions for cosmological statistics such as the power 
spectrum and the mass function eventually requires the treat¬ 
ment of gas physics, star formation, and feedback processes. 
As we explore smaller scales, these effects become more im¬ 
portant. For the matter power spectrum, as an example, we 
expect effects at the 10-20% level at scales of k ~ 10 /iMpc" 1 . 
Accounting for these effects via simulations poses two major 
difficulties: (i) the simulations are computationally very expen¬ 
sive, so only a limited number can be carried out; (ii) the phys¬ 
ical models used in the simulations are too crude to be truly 
predictive; there is a rather large uncertainty in the modeling. 
Steady progress in this area is being made as hydrodynamical 
simulations become more affordable and understanding of the 
influence of baryonic physics improves. Examples of recent 
attempts to mitigate these uncertainties include the alteration 
of the halo concentration-mass relation and exploring how this 
change feeds back in, e.g., the power spectrum (Zentner et al. 
12008) . This approach can be explored within the halo model 
to obtain predictions for the power spectrum. More recently, 
it has been extended to include a range of hydrodynamic sim¬ 
ulations and the calibration of how baryonic physics can bias 
the measurements from surveys like DES and LSST (Zentner 
et al.|20l7}. The idea is not to rely on carrying out a hydrody 


namic simulation that is accurate in all its details, but rather to 
construct an approach that can mitigate uncertainties and biases 
due to baryonic effects in robust ways. 

Finally, a major challenge is posed by the need for accurate 
predictions for covariances. The main difficulty here is the very 
large number of simulations that might be required to deter¬ 
mine covariances at the accuracy needed as well as the uncer¬ 
tainties due to different cosmological inputs. Some early work 
in tackling these questions via emulation has been carried out 
by Schneider et al. ( |2008 1 . Morrison & Schneiderj ( [2017} fo¬ 
cus in particular on the question of cosmology dependence and 
emulation, basing their analysis on a halo model approach. As 
for baryonic effects, more work is required but first results are 
encouraging. We will not address either of these problems here 
but stress that our general emulation approach is directly appli¬ 
cable to these problems. 

To attack the aforementioned challenges a multi-step com¬ 
prehensive program has to be designed and carried out: (i) 
calibrate the underlying dark matter power spectrum and mass 
function for an 8-parameter model space out to k ~ 5Mpc 1 and 
mass between 2 • 1O 12 M0 and 10 I5 M e (at z = 0, smaller cut¬ 
offs for high masses at higher redshifts) at high accuracy; (ii) 
extend the Grange to higher k and mass range to lower masses 
by augmenting the simulation suite with a small number of sim¬ 
ulations with very high mass resolution and use smart extrapo- 
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lation approaches where applicable; (iii) add a small number of 
high-fidelity simulations including gas physics and re-calibrate 
the power spectra and mass functions for these effects; (iv) 
gauge the remaining uncertainties in the predictions and fold 
them into the calibration step. 

Such a program is obviously very ambitious and will require 
significant effort and resources to be completed. In addition, 
community support and input are crucial - the simulation suite 
that will result from such a project can be used for a large num¬ 
ber of scientific projects and has to be made available to the 
entire community. In order to accommodate as many follow-up 
projects as possible, it is important to have most analysis tools 
already in place and a comprehensive plan for post-processing 
so that the relevant data can be properly archived. 

The goal of the current paper is to focus on the first step in 
some detail and demonstrate with surrogate models (such as the 
linear power spectrum and an estimated mass function) and a 
set of high resolution simulations that it is possible to success¬ 
fully complete this program in accordance with the timelines 
of the major surveys. The paper will provide a roadmap that 
delivers prediction schemes at different levels of accuracy at 
different times - culminating finally in high-accuracy emula¬ 
tors that will be required for extracting the most science from 
surveys such as LSST. It introduces a new statistical approach 
for designing a simulation campaign that will lead to improved 
emulators over time. We also show results from simulations to 
investigate the optimal size of the simulations to be carried out, 
considering limited computational resources as well as accu¬ 
racy requirements. We strongly believe that the community has 
to discuss and add to this roadmap to make such a project suc¬ 
cessful. Therefore, a paper outlining the plan and the solutions 
to the problem and showing that it is feasible before embarking 
on this endeavor is crucial. 

The paper is organized as follows. In Section [2j we discuss 
the parameter ranges covered by our study and how we derived 
them. We then introduce a new design strategy for a simula¬ 
tion campaign that can be augmented over time and will lead 
to higher accuracy emulators with each new set of simulations 
(Section [3}. A short overview on the N-body simulations that 
will be employed for the project is presented in Section [4] In 
the following sections, we focus on precision predictions for 
the dark matter power spectrum (Section [5j and the halo mass 
function (Sectional using this new strategy. We present conver¬ 
gence studies and demonstrate that a set of ~ 100 cosmological 
models suffice to reach the required accuracy. We also use a set 
of high resolution simulations to investigate the requirements 
for each individual model to mitigate uncertainties due to real¬ 
ization scatter. In Section [7] we summarize additional emulator 
projects that can be envisioned and conclude in Section [8] with 
a summary and short discussion. 

2. PARAMETER RANGES 


As discussed in the Introduction, a major challenge posed by 
upcoming surveys is the increase in parameter space that must 
be handled. In Heitmann et ali]( |2009j l, we discussed the require¬ 
ments for ongoing weak lensing surveys and concluded that a 
five dimensional parameter space with pg = {u! m ,ujb,ag,n s ,wo} 
was sufficient to capture the information from contemporary 
measurements. In this scenario, we assume flatness, no running 
of the spectral index n s , massless neutrinos, and constraints on 
the Hubble parameter h from measurements of the distance to 
last scattering. 


In this paper, we build upon the results obtained in Heit- 



z 

Fig. 1.— Constraints for the dark energy equation of state parameter w(z ) 
from supernova and CMB data in blue (light blue: 95% confide nce level , dark 
blue: 68% confidence level, red dashed line: w = —1) from |Holsclaw et al.| 
(2011) . The inclusion of baryon acoustic oscillations would tighten the con¬ 
straints further. The red shaded region shows the range for w(z) covered by 
our choices for wq and w a for our main design, comfortably including the best 
current constraints. 


man n et ah] ([2009) and Planck constraints ( |Ade et ah] 12013) 1 
and increase the parameter space by three additional parame¬ 
ters. We extend our previous analysis by letting h be a free 
parameter, allow for a dynamical dark energy component via a 
two-dimensional parameterization (wo, w fl ), and allow for mas¬ 
sive neutrinos via to,, but still keep the assumption of three rel¬ 
ativistic species. For the first five parameters we keep the same 
parameter ranges as in Heitmann et al. [ ( |2009[ ) and we refer the 
reader to that paper for more details on the reasoning behind 
our choices. We provide justifications for the chosen additional 
parameter ranges below. With the Planck data at hand, we keep 
the assumptions of flatness, f \ = 0, and no running of the spec¬ 
tral index, d log n s /d log k = 0. In summary, we investigate the 
following eight parameters within the ranges: 


0.12 <w m < 0.155 

(1) 

0.0215 <co b < 0.0235 

(2) 

0.7 < cr 8 < 0.9 

(3) 

0.55 < li <0.85 

(4) 

0.85 < n s < 1.05 

(5) 

-1.3 < wo < -0.7 

(6) 

-1.5 < w a < 1.15 

(7) 

0.0 < <0.01. 

(8) 


For these parameter ranges our aim is to obtain high preci¬ 
sion predictions for different observables, at the percent level 
of accuracy. 

In Heitmann et al. ( |2009| ), we determined the best fit value 
for the Hubble parameter h for each model from CMB mea¬ 
surements of the distance to last scattering. The power spec¬ 
trum emulator (Lawrence et al.|2010j > simply output the best fit 
value for any other emulated model. In the extended Coyote 
emulator (Heitmann et al. 2014) we provide the power spec¬ 
trum prediction out to higher k- and “-values as well as treat 
the Hubble parameter as a free parameter in an approximate 
way. This extended emulator is based on the same cosmologi¬ 
cal models as the original emulator, which did not allow us to 
provide sub-percent accurate predictions including h as a free 
parameter. Instead, we used results from higher order perturba- 
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tion theory to follow the behavior of the power spectrum under 
changing h (for sufficiently small values of k). 

In the current paper, we keep h as a free parameter embed¬ 
ded in the full design. We choose the range according to the 
minimum and maximum value we found in the model space in 
Heitmann et al. ( 2009| . This range is rather large with respect 
to currently available constraints on h - recent measurements 


Riess et al. (2011)> constrain Hq at the 3% accuracy 


shown in 

level to 73.8 ±2.4 km/s/Mpc, suggesting that measurements at 
the level of 1% accuracy are possible - but we want to ensure 
that the previously investigated parameter space is comfortably 
embedded in the new one. In addition, the rather low value for 
Hq = 67.3 ± 1.2 km/s/Mpc found by the Planck team suggests 
that the error estimate in Riess et al. ( 201 I) might be underesti¬ 
mated. 

Next, we discuss our choices for the range of wq and w a . The 
uncertainty in these parameters, in particular w a , is still very 
large and reducing it is a major target of the upcoming surveys. 
We therefore allow a large range for both parameters, informed 
by constraints from supernova and CMB data. In Fig. [I] we 
show constraints obtained from a non-parametric reconstruc¬ 
tion approach for w(z) using recent supernova and CMB con¬ 
straints (see |Holsclaw et al.||2011| for details). The blue re¬ 
gions show the constraints at the 95% and 68% confidence level 
while red shaded region shows the coverage by our main mod¬ 
els. Current constraints are therefore comfortably covered by 
our allowed parameter ranges. A similar view, but this time 
aimed at current parametric constraints on (iv fl . w a ) is shown 
in Fig. [2] Here we show the coverage of the (wo,w fl ) param¬ 
eter plane and scale the axis approximately following Fig. 36 
in |Ade et al.| ( 2013~| which provides constraints from CMB and 
supernova measurements. Again, with the ranges as chosen, 
we provide good coverage of the currently favored regions of 
parameter space. 

Finally, our choice for the range on the neutrino mass pa¬ 
rameter 'jj y is informed by the constraints available from the 
Planck measurements. In Fig. 3]we show different limits given 
by Planck in combination with other surveys. Note that these 
limits were derived for ACDM models. For our broader model 
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Fig. 2.— Coverage of the (wo, w a ) plane. Red crosses show the design points 
with zero neutrino masses, green crosses the first 26 models, blue stars the next 
set of 29 models, and pink boxes the remaining 46 models. The upper dashed 
line shows the limit wo + w a <0 obeyed by all the models. The area shown 
approximately covers the constraints obtained from Planck data combine d with 
the Union2.1 or SNLS supernova data set (see Fig. 36 in |Ade et al.|2013) . 
Models above the dotted line (wq + w a = 0) are not considered! 
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Fig. 3.— Current upper bounds on from Planck and combinations of 
different cosmological probes for a ACDM cosmology (values from Table 10 
in |Ade et al.|201 3] i. WP stands for WMAP polarization data, the BAO mea¬ 
surements include several different surveys, the 'highL' CMB measurements 
are from the Atacama Cosmology Telescope and the South Pole Telescope. 
The red symbols show the values for in our design for the 101+10 models 
we consider throughout the paper. The first ten models have = 0 to provide 
good coverage of the current Standard Model. The symbols have the same 
meaning as in Fig. [2] 
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space they would be less stringent. We therefore choose to fol¬ 
low the most conservative bound (Planck + WP). The WMAP 
team has considered neutrino mass bounds for wCDM mod¬ 
els and even for those our chosen parameter range comfortably 
covers all currently allowed mass values. 


3. DESIGN STRATEGY 


An important challenge for building efficient emulators in the 
nonlinear regime is overcoming the high cost of accurate simu¬ 
lations. It is not possible to evaluate the simulator over and over 
in a brute force fashion. Given this hard constraint, a judicious 
choice of simulation design (i.e., the simulation suite of mod¬ 
els) is essential. For us, this means varying the eight parameters 
over the ranges specified in Eqns.[T]so that the parameter space 
is optimally sampled. 

The emulation strategy relies heavily on Gaussian process 


(GP) models (see, e.g., Sacks etal. 1989 


Rasmussen & Williams 


2006) to interpolate across the simulation results. Johnson et 


al. (1990[) showed that under certain conditions, designs with 


good space-filling properties (e.g., maximizing the minimum 
distance between the design points - maximin designs) are op¬ 
timal for GPs. This is because a prediction from a GP emulator 
is a weighted average of the data inputs (simulation inputs in 
our case), where the outputs from nearby models have higher 
weights than those far from the prediction location. 

A very important new feature with our simulation design is 
that we aim to progressively fill the model space so that inter¬ 
mediate analyses can be performed while the remaining simu¬ 
lations continue to run. Hence, a simulation design is required 
that simultaneously combines (i) good space-filling properties 
for the chosen number of sampling points and (ii) also provides 
good space-filling properties at intermediate stages of the sam¬ 
pling process. Together (i) and (ii) ensure that the final and in¬ 
termediate stage simulation designs give good coverage of the 
model space so that efficient emulators can be constructed in 
a staged, convergent fashion. To achieve this goal, we use a 
space-filling design based on point lattices. 

It is worth noting that Fatin hypercube designs (McKay et al.) 
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1979[ > are frequently used in simulation based “experiments”, as 
in our previous work. These designs impose one-dimensional 
space-filling. To make them better suited to computer model 
emulation with GPs, additional space-filling criteria, e.g., max- 
imin Latin hypercube designs (Morris & Mitchell 1995]), are 
imposed to improve the design. The nested lattice scheme we 
use here has the added benefit of explicitly enforcing space¬ 
filling properties in intermediate stages of the simulation cam¬ 
paign as well as in the final design. 


3.1. Space-filling lattice designs 

A point lattice A is an infinite set of points in K £/ constructed 
as integer multiples of a set of basis vectors that are given as 
columns of a non-singular generating matrix G e 

A(G) = {Gk : k G Z d } c W 1 , (9) 

writing A(G) instead of A only when the generating matrix is 
needed to differentiate between lattices. The familiar Cartesian 
lattice, for example, has the of-dimensional identity matrix as its 
generating matrix. 

Since a lattice A is a linear transformation of the integers Z d , 
A inherits the integer’s Abelian group structure, i.e., for any 
point Awe have x +A = A. A Voronoi cell is the neigh¬ 
bourhood of points in W 1 that are closer to a particular point 
x G A than to any other point in A. Due to the translation in¬ 
variance of A, its Voronoi cells look exactly the same for any 
x G A. Hence, a lattice basis G that is chosen to give a large 
packing sphere diameter among nearest neighbours (maximin 
distance) or small covering diameter (minimax distance) im¬ 
mediately imposes useful space-filling properties for the entire 
design region Ai as illustrated in Figure [4] ( |Conway & Sloane| 
1999] |Johnson et al.|[T990] l. It is traditional to scale all input 
parameters to the unit interval [0,1] to put all parameters on the 
same footing in the analysis. Thus, it is convenient to write the 
input region as the J-dimcnsional unit cube Ai = [0, \] d with 
volume vol(Ad) = 1. 

A lattice design is the intersection of the lattice and the 
design region (or input region) of interest, possibly shifted 
to achieve useful statistical properties. More formally, X = 
Ai D A + p for input region At is obtained by restricting the 
infinite lattice to the region with a possible shift of A by some 
p G R' / . A straightforward approach for finding a lattice design 
is given in Algorithm[l] 


3.2. Searching for a good lattice design 

The choice of G and p in Algorithm [T] leaves only d 2 + d de¬ 
grees of freedom in optimizing further useful experimental de¬ 
sign criteria. This is important to finding good designs for large 
dimensionality d and simulation experiment runsizes n = A|, 
because without the lattice constraint this would be a d ■ n- 
dimensional, non-convex optimization problem. Furthermore, 
known constructions for good generating matrices G are read¬ 
ily available (Conway & Sloan e~p~999) l, leaving only a suitable 
rotation Q G SO (d) to produce alternative bases QG. 

When a lattice is shifted or rotated, design points can move 
in and out of the fixed design region Ai. The expected run- 
size, h, of a randomly rotated and shifted lattice within the 
unit cube Ad is n(G) = E [lAfAf.QG.p)! : peR‘ i ,Q r Q = l] = 
1 / |det G|. To achieve a certain desired runsize N, at least in ex¬ 
pectation, we scale G such that detG = 1 /N, before performing 
Algorithm [2] 


Algorithm 1 A method to produce a lattice design X(Af, G,p) 
in the region Ai = [0,1 ] d for any generating matrix G and shift 
of origin p. 

1. Determine a finite candidate subset of lattice point in¬ 
dices K' C Z d that produce all design points that are 
potentially contained in [0, l] d . 

For instance, enumerate all integer points in the axis-aligned bound¬ 
ing box that contains the inverse image C' of all corners of the design 
region M, C 1 = { G 1 (c-p): c £ {0, l} rf }. With lower bound 1 £ Z rf 
that is constructed as = inin{/c, : k £ C'} J and analogous upper 
bound u £ Z d , the candidate indices are K' = {/i,Zj + X 

{h,h+ 1, • ■ ■ , 112 } x • • ■ x {lj,ld+ 1,.. ■ ,uj}. 

2. Produce a superset of the design as X' - 
{Gk + p : k G K'}. 

3. Keep all points from X' with coordinates 0 < Xij < 1 to 
obtain the design X. 

Note, to be able to produce lattice designs for larger d it is essential to improve 
the construction of K' in step 1. To sketch the approach here, start by covering 
a one-dimensional projection of G~*(.M —p) to build the design iteratively, 
adding one dimension at a time. 


Algorithm 2 A method to find a shift p and initial rotation Qo 
of a lattice A that produce a design with a desired runsize N = 
1 /detG and other possible design criteria. 

1. Start with an empty set of known designs 1C. 

2. Do the following for a certain number of trials (batch 
size). 

(a) Choose a random shift p G Ai and produce a ran¬ 
dom rotation Q via QR-decomposition of a ran¬ 
dom matrix. 

(b) Generate the lattice design A(A4,QG,p) via Al¬ 
gorithm [T] 

(c) Calculate all design performance criteria c/(X) 
and store them as assessment c. This may just 
include the desired runsize N criterion c: runs i zc = 
Kl-X'l —2V)|. If ||c|| = 0, terminate the search and 
return design X with rotation Qo and p. 

(d) Add the tuple (Q, p, c) to the current batch of trial 
designs. 

3. Add the current batch of designs to the set of known de¬ 
signs 1C. 

Repeat step 2, if the assessments c of known designs 1C 
have changed significantly, for instance, as determined 
by the C\ distance of the percentile functions of em¬ 
pirical distributions of ||c|| in 1C vs. K, joined with the 
current update batch of trials. 

4. Return one of the best known designs, i. e. by choosing 
one with the smallest ||c||. 


3.3. Sequences of nested lattice designs 

A further desirable property of a sequence of design points 
is that they approach a uniform distribution — ideally ensuring 
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good coverage for prediction at earlier run sizes. This means, 
for example, that even with a potential simulation budget of 200 
simulation runs, the first 25 could be chosen so that they can 
inform a decision as to whether resources should be allocated 
to compute a further 25 runs, then 50, and then another 100 
experimental points to gain further insight. 

In this section, we show how to obtain a sequence of lat¬ 
tice designs that progressively fills the parameter (input) space. 
First, observe that a sequence of nested lattices is a hierarchy 
of lattices. To decompose a fine lattice A(G) into coarser lat¬ 
tices, the generating matrix changes into GK for dilation matrix 
K. Central to the construction of nested lattice designs is that 
A(Go) D A(Gi) -£=> Gq'Gi = K £ Z dxd . In other words, any 
sequence of nested lattices 

...dA(G 0 )dA(G 1 )d... (10) 

is equivalent to a sequence of dilation matrices K/ being integer, 
where G/+i = G/K/. The initial lattice A(Go) itself, however, 
does not have to be a subset of the integers. 

The expected runsize n or density of A(G) is increased by a 
factor of |det K| = 3 £ Z as opposed to the coarser sub-lattice 
A(GK). It is possible to choose the factor /I as low as two. So, 
we can construct a sequence of lattice designs with expected 
run-sizes that change by a factor of two. 

Here we provide a new method for constructing a sequence 
of nested lattice designs. This is achieved by developing a con¬ 
struction for a single dilation matrix, K, so that K 1 can be 
used many times to obtain refined designs that yield large run 
sizes n. That is, we start with a coarse lattice, and apply the 
inverted dilation matrix, to get larger designs. The larger lattice 
designs contain the smaller lattice designs and thus have good 
space-filling properties at intermediate run sizes if performed in 
a sequence. 

More precisely, our construction is giving all non-equivalent 
integer matrices K that give a sequence of lattices 

... D A(Gr') D A(G) D A(GK) D ... 


that also have a rotational symmetry 

QG = GK (11) 

for a uniformly a-scaled rotation matrix Q that fulfills Q r Q = 
a 2 I with a d = detQ = detK. 

Equation im suggests an equivalent viewpoint to our con¬ 
struction — the sequence of nested designs is constructed from 
uniformly scaled rotations Q of a specific initial lattice A(G). 
To find a suitable initialization X for the design sequence, the 
generating matrix G can be scaled and rotated to optimize some 
criteria using Algorithm [2] Coarser sub-designs X) based on 
A(Q ; G) or equivalently AfGK 7 ) for / = 1,2,..., as well as re¬ 
fined designs for negative l can be obtained via Algorithm[T] 

A hierarchical ordering of points is obtained by decompos¬ 
ing X and by running the points in X\ without repetition in de¬ 
creasing order starting with the largest /. In practice, the I th 
sub-grid X) is determined by retaining all points of {K 'k : 
k £ G 1 (x - p) : x £ X} that have integer coordinates only. A 
one-dimensional example X\ with K = 2 are the even numbers 
formed from integer numbers that after a multiplication by K 1 
are still integer. 

The choice of generating matrix G and rotation-scale matrix 
Q (Eqn. [JT)i to produce a nested design sequence are not inde¬ 
pendent. Bergner ( |2011} > demonstrated that, given d and sub¬ 
sampling rate 3 = detK, there are a limited number of integer 
matrices K with corresponding rotation-scale matrices, Q, and 



Fig. 4.— The 19 point set (a) is a subset of the 37 point refined design (b), 
which is a rotated version of (a) that has twice the density of sample points. 


G, that can be considered — only four possibilities, K and K r 
with possible scaling by -1, exist in case of odd d and prime 3. 

The lattices obtained in this manner have good space-filling 
properties that lead to efficient prediction for GPs. In addition, 
the sequences from one coarse lattice to increasingly refined 
lattices of the same type gives very predictable improvement 
for the space filling properties (see illustration below). 

Furthermore, unlike the Cartesian lattice where the size of the 
simulation suite grows exponentially with d, the non-rectangular 
lattices we use can achieve virtually any run-size. Perfect ratios 
j3 among runsizes \X/ | are not guaranteed, because some of the 
points in the refined design can lie outside of M. With Al¬ 
gorithm [2] it is possible in step ( |2c| ) to simultaneously match 
runsizes on multiple levels of coarsening l, as additional crite¬ 


ria Cl. 

As an illustration, consider the 2-d example in Figure [4] Op¬ 
timal choices of Q and G are found using the method described 
in Bergner (2011 1 . A coarse design of 19 points is found in 
Figure[4ja). The model runs represented by these points would 
be the first to perform in the simulation and become available 
for intermediate data analysis. The lattice is refined by adding 
another 18 points (the number of simulations is almost doubled, 
but one of the design points in the refined lattice fell outside of 
[0,1 ] d ). This is done by applying the inverse of a dilation ma¬ 
trix to the coarse lattice (alternatively, the dilation matrix could 
be applied to the more dense lattice to obtain the coarse lat- 
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tice). The combined design is shown in Figure |4jb). Notice 
that both the initial lattice and the refined lattice are rotated and 
scaled versions of the same lattice. Furthermore, the shape of 
the Voronoi regions are the same, and the volume of the refined 
cells (area in 2-d) is one-half the volume of the coarse cells. 

The designs that we performed in the 8 -dimensional input re¬ 
gion above were constructed in exactly the same manner. That 
is, an initial lattice design was found for a run size of n\ =25. 
The next ri 2 = 21 points were found by applying the inverse of a 
dilation matrix to the lattice and using the points that intersec¬ 
tion the 8 -dimensional unit hypercube. Finally, the last n 3 = 47 
points were found in exactly the same way. 

To summarize, space-filling design sequences can be ob¬ 
tained by extending a good initial lattice design to a coarsen¬ 
ing or refining sequence of lattices. The rotational symmetry 
of the sequences developed here ensures directionally uniform 
density at each level of refinement. 

3.4. Cosmological Considerations 


from a ACDM simulation suite throughout the following sec¬ 
tions. We give a brief overview of the simulations here. 


4.1. Required Modifications in the N-body Code 

The modifications in the HACC code to include neutrinos 
and dynamical dark energy models are straightforward and 
have been discussed in |Upadhye et al.| ( |2014[ ). In that work, 
HACC simulations were used to test the range of validity of the 
TimeRG (“Renormalization Group”) perturbation theory ap¬ 
proach, introduced by Pietroni (2008) and Lesgourgues et al. 
(2009), focusing on massive neutrinos and dynamical dark en¬ 
ergy models. This work will be an important ingredient for our 
planned power spectrum emulator development. We summarize 
the main points here and remind the reader of the approximate 
treatment of massive neutrinos that we employ. The impact of 
dynamical dark energy and neutrinos are taken into account by 
(i) modifying the initializer and (ii) including both effects in the 
background evolution. 


Currently, observational constraints on neutrino masses are 
only upper bounds, and in the Standard Model of Cosmology 
we still have m v = 0. From a design stand-point, this puts the 
currently best-fit value for one of our parameters at the edge of 
the design (since centering the neutrino mass range around zero 
is excluded). In general, the accuracy of the emulator predic¬ 
tions start to degrade when they approach a design edge, due to 
lack of nearby points. In order to mitigate this problem, we add 
a set of 10 simulations with m v = 0 to the design. We choose the 
values for the remaining 7 parameters according to a symmetric 
latin hypercube design. 

While checking the accuracy properties of a first test emu¬ 
lator, we noticed that models with wo + w„ approaching zero 
were not predicted very well. This was due to the fact that 
the behavior of these models on large scales is different com¬ 
pared to models closer to the Standard Model. With the other 
six parameters fixed, the correlation between two models with 
wo + w„ close to zero is lower than between two models with 
wo +w fl far from zero. Observationally, wo + w a close to zero is 
already disfavored and predictions in this regime do not neces¬ 
sarily need to be accurate at the percent level. Nevertheless we 
aim to provide reliable predictions over the full parameter range 
we cover. Rather than alter our emulation scheme, we introduce 
a transformation of these two parameters to scale the distances 
in these parameters in a way that produces stationary behav¬ 
ior. After some exploratory data analysis, we replaced (wo,w fl ) 
in the design space with (wo,[-wo-w fl ] 1,/4 ). This transforma¬ 
tion results in the desired behavior. This special treatment of 
(wo,w fl ) pairs can be seen in Figure [2] — the design provides 
slightly more points close to the limit wo + w a = 0 . 

4. JV-BODY SIMULATIONS 


The parameter range outlined in Section [2] requires some 
extensions of the standard setup for wCDM N-body simula¬ 
tions. In particular, dynamical dark energy models defined by 
an equation of state parametrized by (wo,w fl ) have to be in¬ 
cluded as well as a treatment of massive neutrinos. In the fol¬ 
lowing we will give a brief description of our implementations 
of these additions in the A'-body code HACC (Hardware/Hybrid 


Accelerated Cosmology Code) described in detail in Habib et 
laL] d 20 l 6 | ). 


In order to estimate the overall computational needs with re¬ 
gard to box size, mass and force resolution, number of real¬ 
izations, and other considerations, we will also discuss results 


4.1.1. Neutrinos 


Simulating neutrinos fully self-consistently as a separate par¬ 
ticle species in a cosmological simulation is nontrivial mainly 
due to two reasons: (i) the very large neutrino velocities at early 
times due to the thermal velocity contribution which make very 
small time steps necessary, (ii) the large mass difference be¬ 
tween the tracer particles representing dark matter and neutri¬ 
nos respectively, which can lead to scattering artifacts if not 
treated properly. Over the years, many different solutions have 
been suggested to overcome these problems, from starting the 
(neutrino) simulation very late (leading to a different set of 
problems and inaccuracies) to simulating the neutrinos on a 
coarser grid to treating the neutrinos only linearly. An in¬ 
complete list of papers employing different solutions include 


Klypin et al. (1993), Gardini et al. (1999), Brandbyge et al. 

(2008), 

Brandbyge & Hannestad ( 

2009), 

Brandbyge & Hannes- 

tad (201 

0 

), Viel et al. (2010), Birc 

et al. ( 

2012), and Agarwal & 

Feldman 

2011 ). 


In our work, we follow closely the approach used in |Agar-| 
|wal & Fel dman (2011) with some modifications. We do not in¬ 
clude the nonlinear evolution of the neutrinos but instead evolve 
the cold dark matter-baryon component (cb) only, including the 
neutrinos in the background evolution. In order to obtain the to¬ 
tal matter power spectrum, we add at any given output the linear 
component of the neutrino power spectrum to the nonlinear cb 
matter power spectrum in the following way: 




yjp^(k,a)+yjp^(k,a) 1 . 


( 12 ) 


This treatment is justified since neutrinos are massive enough 
to behave matter-like close to z = 0 but light enough to not clus¬ 
ter strongly on small scales. In Upadhye et al. (2014 1 we in¬ 
vestigated the validity of this approximation on large to quasi- 
nonlinear scales and found sufficiently small errors for neutrino 
masses below the currently excluded value. 

When setting the initial conditions, we include the neu¬ 
trino contributions in the transfer function and generate a lin¬ 
ear power spectrum at z = 0 normalized to a given a%. Then, 
the cb component of the matter power spectrum is moved back 
to the initial redshift with a scale-independent growth function 
(the growth function that includes neutrinos would be scale- 
dependent). This growth function includes all the species in the 
homogenous background. Note that using a scale-dependent 
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growth function would be inconsistent since the evolution to 
z = 0 is carried out for the cb component. The setup described 
here ensures that at z = 0 the power spectrum matches the full 
linear power spectrum on large scales. Finally, in the evolution 
the massive neutrinos are included in the background evolution 
equation, but we do not source the neutrino growth in the Pois¬ 
son equation. For a detailed discussion of the implementation 


as well as various tests, the reader is referred to Upadhye et al. 
{2014} . 


4.1.2. Dynamical Dark Energy Equation of State 

Next, we discuss our treatment of a dynamical dark energy 
equation of state. Without a compelling theory for a dark en¬ 
ergy model beyond the cosmological constant, the simplest way 
to express a dynamical origin for dark energy is to parameterize 
the dark energy equation of state w(z) via a time varying func¬ 
tional form described by two parameters (wo, w a ). A commonly 
used form is that suggested by|Chevalier & Polarski|(|2001J and 
|Linder| ( |2003] >: 

w(a) = WQ + w a (l-a). (13) 

In order to account for the time dependent dark energy equa¬ 
tion of state in the initial conditions of the simulations, we need 
to generate a transfer function that includes the effects from 
w(z). For this, we modify the publicly available version of 
CAMB ( jFewis et al.{2000) and include the dynamical dark en¬ 
ergy in the evolution of the background as well as modify the 
equations describing the density and velocity perturbations in 
the dark energ) 

To modify the background cosmology, we simply have to 
change the equation describing the evolution of conformal time 
as a function of the scale factor. Modifying the equations de¬ 
scribing the perturbations of dark energy requires an expres¬ 
sion for the speed of sound. In analogy with a simple scalar 
field model, we assume that this is the speed of light. As a re¬ 
sult, perturbations of dark energy develop only on the Hubble 
scale. A second issue is that the equations describing the evo¬ 
lution of the velocity perturbations include terms of the form 
c 2 /( 1 + w(a)), where c s is the speed of sound in the rest frame 
of dark energy. For those values of (wo, w„) for which the equa¬ 
tion of state passes through -1, this results in a singularity. We 
treat this by assuming the microscopic properties of dark en¬ 
ergy do not modify the power spectra, and so this term is small 
when w{a) crosses -1. Next, we need to include (wo,w fl ) in 
the growth function in order to set up the initial conditions at 
the desired redshift. This is achieved by evaluating the inte¬ 
gral over w(a) that appears in H 2 (a ) in the growth function: 
-3 /, “da'[w(a')/a r ] = exp[-3w a (l - a)]. Finally, the 

modifications in the Vlasov-Poisson equation for the time evo¬ 
lution to include (wo,w fl ) are also simple. In the case of HACC, 
the expansion factor a is used as the time variable, and there¬ 
fore the only change again is in I Hz.) and the corresponding 
scale factor expressions. 


4.2. A CDM Simulation 

In order to test our overall simulation strategy with regard 
to resolution and sufficient statistics, we carry out two high- 
resolution ACDM simulations. These simulations represent the 
specification of the main simulation suite that will be carried 
out for each model in our future work. In addition to these high 

3 The modified CAMB version can be downloaded at: 
http://www.hep.anl.gov/cosmology/pert.html 


Table 1 

Simulation specifications 


n p 

L[Mpc] 

Realizations 

Force resolution [kpc] 

3200 3 

2100 

1 

6.6 (high-res) 

4096 3 

5634 

1 

13.8 (high-res) 

512 3 

1300 

16 

1270 (PM) 


mass and force resolution simulations, we will augment each 
model with a set of lower resolution simulations to cover the 
full dynamic range with sufficient accuracy as required. The 
cosmological parameters are chosen to be close to the best-fit 
WMAP-7 parameters (Komatsu et a l,|2011 i: uj m = 0.1335, w* = 
0.02258, n s = 0.963 , w= —1.0, erg = 0.8 ,h = 0.71. In one simu¬ 
lation we evolve 3200 3 particles in a (2100 Mpc ) 3 volume, lead¬ 
ing to a mass resolution of m p ~ 10 10 Mq. The force resolution 
is chosen to be ~ 6.6 kpc. This simulation has been previously 
used to validate the accuracy of a wCDM power spectrum em¬ 
ulator (Heitmann et a l,|2014[) and to create an emulator for the 
galaxy power spectrum ( Kwan et al.|2013a| > based on halo oc¬ 
cupation distribution (HOD) modeling. The idea behind using 
an HOD to generate galaxy catalogs is straightforward (Kauff- 


mann et al. 1997||Jing et al. 

1998 ; Benson et al. 2000; Peacock 

& Smith 2000. Seljak 2000 

Berlind & Weinberg 2002; Zheng 


population (central and satellites) for a halo as a function of 
halo mass. The HOD parameters are obtained from matching 
the resulting galaxy correlation function to observations. For 
different types of galaxies, the HOD will be specified by differ¬ 
ent parameter values. The resolution of our simulations is suf¬ 
ficient for a wide range of investigations and to build detailed 
synthetic skies in different wavebands. 

The specifications of the additional simulations needed to 
generate accurate emulators depend on the statistics of inter¬ 
est. As we show below, for the power spectrum predictions 
we require a set of particle-mesh (PM) runs for the quasi-linear 
regime to suppress the noise due to realization scatter. For the 
mass function, we require large volume simulations with high 
force resolution in addition to the simulation described above 
to cover the cluster mass range. As we show in Section [ 6 j a 
simulation covering a volume of (5600 Mpc ) 3 evolving 4096 3 
particles is optimal for this purpose. Due to the modest mass 
resolution of ~ 1 O U M 0 and - as we show below - the rather 
small number of time steps needed to achieve sufficient accu¬ 
racy, the cost of these large simulations is reasonable. A sum¬ 
mary of the ACDM simulations used in this paper is given in 
Table □ 


5. THE POWER SPECTRUM 


5.1. Smooth Power Spectrum Predictions 


Following the discussions in Lawrence et al. (2010 1 and |Heit-] 
man n et al.| ( |2014| l, a major obstacle in extracting a smooth 
power spectrum from a simulation campaign is due to realiza¬ 
tion scatter. In order to solve this problem, we have shown that a 
smoothed power spectrum can be achieved by matching higher 
order perturbation theory on the largest scales with medium 
resolution, particle mesh (PM) runs (many realizations) on in¬ 
termediate scales, and a high resolution simulation on small 
scales. In addition to this matching procedure, we have de¬ 
veloped a sophisticated smoothing method based on a process 
convolution approach described in Higdon (|2002[l. The main 
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Fig. 5.— Upper panel: Power spectrum at z = 0 assembled from RPT (green 
dashed), the average of 16 PM simulations (red dashed), and one high reso¬ 
lution simulation at large scales (light blue dotted). The black curve under¬ 
neath the simulation/perturbation theory results shows the smooth power spec¬ 
trum extracted from this ensemble. In order to demonstrate the quality of the 
smoothing procedure, we show the ratio of these results in the lower panel, 
including error bars calculated from the finite number of modes at each bin for 
the simulation results. 


idea behind this approach is to build a smooth function as a 
moving average of a simple stochastic process. The moving 
average uses a smoothing kernel with a varying width - in re¬ 
gions with distinct features like the baryon acoustic oscillations 
the window is narrow while in very smooth regions such the 
very small scales, the window can be wider. We will use the 
same method again here. We cover a k range between 0.001 to 
6 Mpc 1 via matching (i) renormalized perturbation theory out 
to k ~ 0.04 Mpc -1 for z = 0, 1 and k ~ 0.14 Mpc -1 for z = 2; 
(ii) 16 realizations of PM simulations with 512 3 particles on a 
1024 3 grid in a (1300 Mpc) 3 volume out to 0.25 Mpc -1 ; (iii) 
a (2100Mpc) 3 volume, high-resolution simulation with 3200 3 
particles and 6.6 kpc force resolution to cover the remaining 
k-range. 

Figure [5] demonstrates the matching and smoothing strategy 
for our fiducial ACDM model, M000, described in Section [P] 
The power spectrum was obtained from the simulation in the 
same way as described in Heitmann et al. (2010) via mapping 
the particles onto a large grid (6400" in this case) using a cloud- 
in-cell (CIC) assignment, applying an FFT, correcting for the 
charge-assignment window function and averaging the results 
in bins in jk|. We show results for both P(k) and the dimension¬ 
less power spectmm A 2 (k), connected via the following rela¬ 
tion: 

k 3 P(k ) 


A z (k): 


27r 2 


(14) 


We only present results at z = 0 here, for z, = 1 and z = 2 the 
reader is referred to the Appendix of Heitmann et al.| (j2014| . 
The upper panel shows the full power spectrum assembled from 
RPT, PM simulations, and one high-resolution run on top of the 
smoothed power spectrum. On the logarithmic scales used here, 
even without the smoothing, the result is already very good. 



Fig. 6.— Power spectrum model space covered by 101 cosmological models 
within the range described in Eqs. jT|. shown in red. Test models are shown in 
blue, where the light blue curve shows M000, a ACDM model, and the dark 
blue curves show the additional models specified in Table[2] 


The lower panel shows the ratio of the smoothed power spec¬ 
trum to the simulations; the aim of this figure is to demonstrate 
that the matching does not lead to any bias. These results show 
that our matching and smoothing strategy works as desired. We 
also point out that the large simulation volume paired with good 
mass resolution allows us to push out to higher k without hav¬ 
ing to rely on small simulation volumes, which can induce some 
inaccuracy, as was shown in Heitmann et al. ( 2014j >. 


5.2. Emulator Performance 

Following the strategy outlined in Section [3] we now discuss 
the accuracy we expect to achieve from 101 simulation models, 
obtained in a three-step simulation campaign. The first set of 
power spectra covers 26 models, the second set 55, and the final 
set covers all 101 models. In addition, as discussed in Section[3] 
we add 10 models with in,, = 0 to properly cover this important 
edge of the design hypercube. 

For our tests, we use linear power spectrum predictions gen¬ 
erated with CAMB. This allows us to obtain many power spectra 
with nominal computational effort. At the same time, the linear 
power spectrum has been proven to be a very good testbed for 
the final nonlinear emulator (see, e.g., Heitmann et ak] |2009) l. 
Figure [6] shows the model space we cover via the cosmological 
parameter ranges given in Eqs. G and the cosmological models 
used to test the accuracy of the emulator. The model specifica¬ 
tions are given in Table[2]in Appendixes models E001 to E010. 
These extra models were chosen by generating another design 
within the parameter hypercube. We do not use a model on the 
edge of the hypercube for our test. At the edges, the accuracy 
requirements are not as stringent as those models are obser- 
vationally already ruled out. Nevertheless, for the purpose of 
an inverse analysis using Markov chain Monte Carlo (MCMC) 
methods, predictions over a wide parameter space are desirable. 

Results from the accuracy test are displayed in Figure [7] We 
show the error of the emulators obtained from 26, 55, and 101 
models in the left three panels and the results including the ex¬ 
tra 10 runs for m v = 0 in the three panels on the right. The 
accuracy improves considerably on adding more cosmological 
models particularly on small scales. Considering the large pa¬ 
rameter space we cover, the accuracy obtained from only 26 
models is already extremely good. The two models that are not 
predicted well on large scales (out of 26) are E007 and E009. 
For E007, we have wq + w a = -0.011 and for E009 we have 
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Fig. 7.— Emulator test: An emulator for the linear power spectrum is tested 
using 10 reference cosmological models. The emulator predictions built from 
up to 101 models are compared to the exact linear spectra from CAMB in the 
left panels and predictions from up to 111 models in the right panels. Shown 
is the ratio of prediction to “truth”, the emulator predictions are accurate at the 
1% level (dashed line) on small scales. The upper panels shows the prediction 
for an emulator based on 26/36 models, the middle panel is based on 55/65 
models, and for the predictions shown in the lower panel all 101/111 models 
were used. The red line shows the prediction accuracy for a ACDM model. 


wo + w a = -0.0222. Both of these models therefore fall in the 
class of so-called “early dark energy models”. The reason for 
the inaccurate prediction is a slight increase in power on large 
scales in these models. Nevertheless, after increasing the model 
space to 56 and 101 the results of these models are also accu¬ 
rately captured. Finally, in addition to the ten models we also 
show the prediction for the ACDM case in red. Here one can 
see how the additional 10 models with zero neutrino mass help 
to improve the prediction accuracy. Overall, with the complete 
set of models, we are able to provide predictions at the ±1% 
level (shown by the dashed lines in all panels) over all 8 pa¬ 
rameters out to k ~ IOMpc 1 . Extending the range to higher 
k values involves running higher-resolution simulations (both 
force and mass) and properly accounting for baryonic effects. 
As discussed in the Introduction, results from initial investi¬ 
gations are becoming available for the latter. In the future, a 
combination of extrapolation to large k values and modeling of 
baryonic effects will be required for robust predictions on the 
smaller length scales relevant to cosmological analyses. 

6. THE MASS FUNCTION 
6.1. Smooth Mass Function Prediction 

In the same way as for the power spectrum, building a mass 
function emulator requires a smooth prediction for the mass 
function for each cosmological model in our simulation de¬ 
sign. For the mass function, this task is much easier than for 
the power spectrum, since no subtle features such as baryonic 



Fig. 8.— Upper panel: Mass function including Poisson error bars from 
the Mira-Titan Universe simulation (red) to cover the low mass end of the 
mass function and an additional large volume simulation, spanning a (5664 
Mpc) 3 volume (blue) to cover the clu ster reg ime. Lower panel : Ratio of the 
simulation results with respect to the |Bhattacharya et al. 1(2011) fit, assuming 
universality. The yellow shaded region shows the 2% error band. 


wiggles are present. The much more challenging aspect here 
is to obtain enough statistics with regard to halo counts; as has 
been shown previously (see, e.g., |Lukic et al.||2007[ |Tinker et| 


al.|2008 Crocce et al.|2010] Bhattach arya et al.|201 ![ > captur- 

ing the mass function accurately at the cluster mass scale is not 
easy. This is due to the exponential drop-off of the mass func¬ 
tion at large mass. Small uncertainties in the mass estimates for 
massive halos are amplified and lead to effects at the level of 
tens of percent in the mass function itself. In order to generate 
an accurate mass function from simulations, it is important to 
have (i) a large enough volume to capture a large number of 
cluster sized halos, (ii) good force resolution to ensure that the 
halo masses are accurate, and (iii) good mass resolution to en¬ 
sure that halos are sufficiently resolved to ensure accurate mass 
determinations. 

The mass range of interest for cosmological surveys in the 
area of “cluster cosmology” actually covers the range from 
galaxy groups to large clusters. With a mass resolution of 
m p ~ 1.05 • 10 10 Mg for the ACDM case, the Mira-Titan Uni¬ 
verse comfortably resolves this range of masses. With a vol¬ 
ume of (2100Mpc) 3 the number of very massive clusters is not 
very large, and we therefore limit the use of the Mira-Titan 
Universe run to a mass range between mhaio ~ 5 • 10 12 Mg to 
wihaio ~ 1 • 10 14 Mq. For the small mass end, halos are there¬ 
fore resolved with ~ 500 particles which minimizes the need 
for (FOF) mass corrections due to undersampling of the halos. 
For the high-mass halos, this cut ensures that we have tens of 
thousands of halos in the relevant mass bins. In order to cover 
the upper end of the mass function, we employ another set of 
simulations. These simulations cover volumes of ~ 5600Mpc 3 
and evolve 4096 3 particles, leading to a mass resolution of 
m p ~ 10 n Mg. These simulations were designed to generate 
mock catalogs for surveys such as BOSS and DESI (Sunayama 
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Fig. 9.— Upper panel: Combined sim ulation result (red) for /(<r) compared 
to the best fit from Bhattacharya et al. (2011) and the best fit to the simula¬ 
tion result. The two fits are basically indistinguishable. Lower panel: Ratio 
between the two fits. The agreement is better than one percent. 


et al. in preparation). At masses of Whaio ~ 1 ■ 10 14 Mq halos are 
resolved with ~ 1000 particles. If we restrict the mass range to 
ffihaio ~ 10 I5 Mq, the last mass bin has more 10,000 halos, lead¬ 
ing to small statistical errors. 

Figure [8] shows the results for the mass function from the 
Mira-Titan Universe simulation for the lower mass end and the 
results from the large simulation for the cluster regime. We use 
the friends-of-friends (FOF) definition for determining the halo 
masses using a linking length of b = 0.2 (Davis et al.| [1985| >. 
This mass function definition has been studied very thoroughly 
in the past and allows us to compare our results to other work. 
(Our discussion can easily be extended to overdensity masses 
as well.) The lower panel displays the ratio of the simulation 
results with respect to a mass function fit developed in Bhat- 
|tacharya et al.[ ( {201 1) . The simulations agree with this fit to 
better than 2% (yellow band) over the full mass range. 

The functional form of the Bhattacharya et al. (2011 1 fit is 
inspired by the original Sheth-Tormen mass function fit ( Sheth 
|& To rmen 1999) but has one additional parameter and allows 
for explicit redshift dependence: 


/ Bhatt (cr,z) =A\j — exp 


2a 2 


1 + 


a%) 


P 1 


$cVa 


05) 

with the following parameters (note that all parameters are dif¬ 
ferent from the original Sheth-Tormen choices): 

0 TTT 0 78R 

A = d^or ; fl= (T^aor ; P = 0 - m ’ « = 1 * 795 ’ (16) 

with 6 C = 1.686. Figure [9] shows the results for f(a(M)) as a 
function of 1 / a(M ) as measured by the combined simulations 
and the fit from Bhattacharya et al. ( 201 lj l. We have refitted 
all parameters for z = 0 using our current simulations and find 



le+14 le+15 

M[M S ] 

Fig. 10.— Mass function error in the regime of cluster masses due to coars¬ 
ening the time stepping. Shown is the ratio of a simulation with 2 sub-cycles 
per long-range step over 5 sub-cycles. In the mass regime where the large vol¬ 
ume runs will be used, the difference is at the sub-percent level. Given the 
computational saving of a factor of 2.5, this small inaccuracy is tolerable. 


extremely good agreement with the previous results: 


A = 0.3315±0.0003; a = 0.7895±0.0085; (17) 

p = 0.8148 ±0.0451; q = 1.8001 ±0.0425, 


shown in Figure [9] The agreement between the fits at the sub¬ 
percent (the lower panel in Figure [9] shows the ratio of the two 
fits) is rather remarkable and reassuring, given the fact that the 
simulations were carried out with different codes. Since the 
simulation results are very smooth already, we do not need to 
employ a special smoothing procedure. Instead, we will be able 
to build the emulator directly from the binned mass function 
datasets, eliminating one source of uncertainty in the emulator 
construction. 

Next, we briefly discuss how to reduce the computational 
costs for the large simulation volume runs. In Sunayama et al. 
(in preparation) a detailed study is carried out concerning the 
number of time steps needed to obtain accurate predictions for 
halo statistics. In the HACC code the time stepper is divided 
into long time steps for the long-range force solver and each 
long time step is subdivided into shorter time steps with short- 
range force solves in order to resolve small scale structures ac¬ 
curately. At a mass resolution of m p ~ 10 u Mq and with the 
mass function as the observational target, the number of sub¬ 
cycles can be reduced without losing accuracy. Figure|T0]shows 
the ratio of the mass function in the mass regime of interest for 
the large simulation with 2 and 5 sub-cycles (earlier tests shown 
in, e.g., Habib et al. ( 2016| > have demonstrated convergence for 
5 sub-cycles). The difference is at the sub-percent level, with 
computational costs reduced by a factor of 2.5. At the same 
time, these simulations will be still very useful for generating 
synthetic sky catalogs for large scale stmcture surveys such as 
DESI. 


Our final mass function emulator will go beyond a single 
mass definition. We will measure mass functions from our sim¬ 
ulations for both FOF and overdensity definitions, varying link¬ 
ing length and overdensity specification. This will make it more 
convenient to compare the results to observational mass mea¬ 
surements and proxies as well as to study connections between 
the different halo definitions. The extension of our results 
shown above to build a flexible emulator like this is straight¬ 
forward. 

Finally, we comment on the higher cluster mass range of the 
mass function. In this paper, we have been very conservative 
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Fig. 11.— Coverage of mass function model space with the simulation 
design described in this paper at z = 0 (red). The dark blue lines show the test 
models and the light blue line the fiducial ACDM mass function. All models 
are based on assuming universality is valid over the full parameter range. This 
assumption is valid at the ~ 10% level of accuracy. 


with regard to pushing the mass range too far. In future work, 
following the discussion in |Lukic et al.| ( [2007| , we will include 
careful investigations on finite volume effects in our simula¬ 
tions allowing us to push to higher masses. We will provide a 
more detailed discussion on this topic in a forthcoming paper, 
where we focus solely on the mass function. 


6.2. Emulator Performance 

In order to test the emulator performance following the de¬ 
sign strategy outlined in Section [3j we need a mass function 
surrogate - just as we used linear theory for the power spec¬ 
trum. First suggested in Jenkins et al. ( |2Q0l) >, writing the mass 
function in terms of o(M) allows for an almost universal de¬ 
scription of the friends-of-friends (FOF) mass function for a 
linking length of b = 0.2 that depends only on the linear power 
spectmm. Over time this universality has been investigated in 
detail (see, e.g., Reed et al.|2007[ Cohn & White|2008) Lukic 


etakj [2007 |Tinker et al.||2008HCrocce et al.||2010UCourtin et 


al.|2010 Bhattacharya et al.|201 1 1 with the conclusion that for 

the specific definition of the b = 0.2-FOF mass function, uni¬ 
versality holds at the 10% level accuracy over a wide range of 
cosmologies. We use this result in the following to investigate 
the accuracy of a mass function emulator over the parameter 
space discussed in this paper. 

We use the universal form of the mass function given in 
Eq. ( fi~6| by Bhattacharya et al. ( |2011[ > to generate mass function 
predictions for all our models. Figure[TT|shows the model space 
covered by the simulation campaign described in this paper as 
well as the test cosmological models used to gauge the emula¬ 
tor accuracy. The ACDM model used throughout the paper is 
also shown. 

Based on these mass function predictions, we build an em¬ 
ulator in the same way as for the power spectrum. Figure [12] 
shows the results for the emulator predictions for the test mod¬ 
els at z = 0. The first two panels show the results from an emu¬ 
lator based on 26 and 36 models. The accuracy of the emulator 
is already within 5%, a very encouraging result. Increasing the 
number of models to 55 (65 including the m v = 0 models) re¬ 
duces the prediction error to ~ 3%, and the final results for 101 
and 111 models promise to yield inaccuracies at the 1% level 
over an 8-dimensional parameter space. This leads us to the 
conclusion that if the mass function prediction from the sim- 




Fig. 12.— Emulator test: An emulator for the mass function is built from 
up to 101 and 111 models, assuming universality. Next, the mass functions 
for 10 additional models are predicted with the emulator and compared to the 
prediction assuming universality. Shown is the ratio of prediction to “truth”, 
the emulator predictions are accurate at the 1% level (dashed line) over the 
full mass range once 101 (111) models are considered. The upper panels show 
the prediction for emulator based on 26 and 36 models, the middle panels are 
based on 55 and 65 models, and for the predictions shown in the lower panels 
all 101 and 111 models were used. As before, the right panels include the 
addition 10 simulations with m v = 0. 


ulation can be provided at sufficient accuracy, building highly 
accurate emulators will be straightforward. 

7. BEYOND THE POWER SPECTRUM AND MASS FUNCTION 

7.1. Emulators 


In the previous sections we have shown examples of two pre¬ 
cision emulators to be built from the Mira-Titan Universe cam¬ 
paign targeted to impacting analysis of data from ongoing and 
upcoming dark energy surveys. The outlined simulation cam¬ 
paign lends itself to building many more emulators. One such 
example is a concentration-mass ( c-M ) emulator. In Kwan et 
|ak| ( |2013a) > it was shown that an accurate c-M emulator can be 
built from a small set of simulations - the paper was based on 
an augmented set of Coyote Universe simulations. The Mira- 
Titan Universe volume and mass and force resolution is suffi¬ 
cient to generate predictions for the c—M relation over a mass 
range covering groups to cluster scales. Figure [13] shows the 
c-M relation measured from the Mira-Titan Universe run at 
z = 0. Given the intrinsic scatter in the c-M relation, the statis¬ 
tics obtained from the Mira-Titan Universe simulation are more 
than sufficient to derive accurate c-M relation predictions over 
the mass range shown. Each halo here is resolved with at least 
1000 particles to provide well-sampled halos when determining 
the concentration (for details about the concentration measure¬ 
ments, see [Bhattacharya et al.|2013} . For comparison, we show 
the best fit power law derived from the Q Continuum simu- 
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Fig. 13.— Concentration-mass relation as measured from the Mira-Titan 
Universe simulation (red points). The red line shows the best-fit power law 
through the simulation points. In bl ue we show t he best fit c—M relation found 
in the Q Continuum simulation {Heitmann et al.|2014| , a simulation with al¬ 
most 100 times better mass resolution, allowing the derivation of a fit out to 
much smaller halo masses. The yellow shaded region shows the intrinsic scat¬ 
ter. Within this scatter, the predictions from the Mira-Titan Universe simula¬ 
tion agree very well with the high-mass resolution Q Continuum simulation. 


lation Heitmann e t alT] ( j2014{ ), a simulation with much higher 
mass resolution. The agreement is very good. 

Another example is an emulator constructed for the galaxy 
power spectrum. In Kwan et al. ( 2013b| ), the Mira-Titan Uni¬ 
verse simulation was used to build a galaxy power spectrum 
emulator covering a wide range of HOD models. With the com¬ 
plete suite of models we will be able to extend that work to 
not only cover different HOD models but to also provide pre¬ 
dictions across cosmologies. Extending our work into redshift 
space is another obvious direction, as redshift space distortion 
power spectrum and correlation function emulators can be eas¬ 
ily obtained based on the simulations outlined here. 


7.2. Sky Maps 

Beyond emulators, the simulation suite described here will 
also be very valuable for producing synthetic sky maps. For 
the LSST Dark Energy Science Collaboration we have gener¬ 
ated an HOD-based galaxy catalog for the large scale structure 
working group. The size and resolution of the main simulations 
will also be sufficient to build weak lensing maps for surveys 
such as DES. Maps tracing the kinematic Sunyaev-Zel’dovich 
effect are currently being developed. We are planning to pub¬ 
licly release these maps to the community. 

8. CONCLUSION 

In this paper we have introduced the Mira-Titan Universe, a 
new simulation campaign to provide accurate predictions for a 
range of cosmological observables and for generating sky maps 
in different wavebands. This work is a major new extension of 
our previous Coyote Universe project in several key respects. 

1. The cosmological parameter space covered is extended 
by three parameters, w„, h, and m v , all of which will 
play a central role in future analyses of dark energy sur¬ 
veys. We have shown that despite the extended parame¬ 
ter space the number of models needed to build reliable 
emulators is still manageable. 

2. A new design strategy has been implemented using 
space-filling lattices. This strategy allows us to improve 


the emulator quality by systematically adding more sim¬ 
ulations to an existing latice design. Our new strat¬ 
egy has the benefit of explicitly enforcing space-filling 
properties in intermediate stages as well as in the final 
design. We would like to stress that this procedure is 
highly nontrivial - extensive tests showed that simply 
adding independent space-filling designs (as in conven¬ 
tional Latin hypercube sampling) will not lead to the 
desired results. 

3. We have extended our investigations beyond the power 
spectrum and demonstrated that our design will also 
provide high-accuracy emulators for the mass function. 
In addition, the simulation set will allow us to build a 
range of other emulators, for the halo concentration- 
mass relation, galaxy correlation functions and power 
spectra, and redshift space distortions. 

4. We have demonstrated with a set of high resolution sim¬ 
ulations that the needed requirements for such a simu¬ 
lation suite can be met. These simulations can then be 
also used for a range of other investigations, in particu¬ 
lar building mock catalogs for different wavebands. 
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APPENDIX 

In this Appendix we list the test cosmologies that were used 
to verify the accuracy of the emulators discussed in the paper. 
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